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This paper deals with the accurate microscopic simulation and macroscopic modeling of 
extreme non-equilibrium phenomena, such as encountered during hypersonic entry into a 
planetary atmosphere. The state-to-state microscopic equations involving internal excitation, 
de-excitation, dissociation, and recombination of nitrogen molecules due to collisions with 
nitrogen atoms are solved time-accurately. Strategies to increase the numerical efficiency are 
discussed. The problem is then modeled using a few macroscopic variables. The model is 
based on reconstructions of the state distribution function using the maximum entropy 
principle. The internal energy space is subdivided into multiple groups in order to better 
describe the non-equilibrium gases. The method of weighted residuals is applied to the 
microscopic equations to obtain macroscopic moment equations and rate coefficients. The 
modeling is completely physics-based, and its accuracy depends only on the assumed 
expression of the state distribution function and the number of groups used. The model 
makes no assumption at the microscopic level, and all possible collisional and radiative 
processes are allowed. The model is applicable to both atoms and molecules and their ions. 
Several limiting cases are presented to show that the model recovers the classical two- 
temperature models if all states are in one group and the model reduces to the microscopic 
equations if each group contains only one state. Numerical examples and model validations 
are carried out for both the uniform and linear distributions. Results show that the original 
over nine thousand microscopic equations can be reduced to 2 macroscopic equations using 
1 to 5 groups with excellent agreement. The computer time is decreased from 18 hours to less 
than 1 second. 


I. Introduction 

D URING hypersonic entry into the atmosphere of a planet or moon, space vehicles may encounter high 
temperature and/or rarified gas environments, such that the flows around them could be significantly affected 
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by various non-equilibrium phenomena. On a macroscopic level, these may be manifested as diffusion, viscous 
stress, heat conduction, internal energy redistribution, chemical reactions, or thermal radiation. They are due to 
excitation of translational and internal modes of gas particles resulting in changes of the microscopic velocity and 
quantum energy level distributions from their equilibrium values. The exact treatment of these phenomena can be 
obtained at the microscopic level by considering the state-to-state collisions between atoms, molecules, and 
electrons, and the emission and absorptions of photons. The microscopic state distribution functions of gas particles 
are governed by generalized Boltzmann equations [1] encompassing the physical, velocity, and internal energy 
spaces. Solving these equations is computationally impractical at the present time even with the most advanced 
supercomputers, since they involve a huge number of quantum energy states and an extremely large range of 
velocities of individual particles. These microscopic equations are therefore modeled by macroscopic equations 
involving a small number of macroscopic unknowns. The accurate modeling and simulation of these non- 
equilibrium processes was started 100 years ago [2, 3], and today it still remains one of the greatest challenges in gas 
dynamics. 

Before the advent of today’s vast computing power, all modeling relied on analytic methods. This required 
many simplifications to the true physics. The main simplification was that the translational and various internal 
modes were all uncoupled. The first internal mode to be modeled was vibrational non-equilibrium due to collisional 
processes by Landau and Teller [4]. It is based on a number of assumptions, which are listed below: 

1) Vibration and rotation are assumed uncoupled and the latter is in translational equilibrium. 

2) Vibrational energy levels are those of a linear harmonic oscillator. 

3) Transitions occur only between immediately adjacent states. 

4) Forward and backward rates are related by the principle of detailed balancing. 

5) Rate constants are proportional to the quantum number of the upper state. 


d ^ ^ ^ 

The model has the form of linear relaxation — - = — . This remarkably simple equation shows that the 

dt x 

vibrational energy e v of the system of oscillators will always tend toward the equilibrium value e * v . However, when 
applying to realistic molecular gases, we observe several shortcomings: 


1) The uncoupled assumption would result in an inaccurate calculation of the internal energy, and that may 
lead the system to relax to a wrong equilibrium value e * v . Reference 5 showed the error in total internal 
energy to be approximately 4.3% at 10,000 K and 11% at 20,000 K compared with accurate quantum 
mechanical calculations. Conversely, errors in temperature for given e* v could be even larger. 

2) In deriving the equation no assumption was ever made as to how the vibrational energy e v is related to a 
temperature. In fact, the vibrational energy was defined precisely as the sum of the vibrational energy of 
all particles, e v = 2 n V ' , where n! is the population of the energy state e l . The population n l can be 

obtained easily by solving the master equation involving only a tridiagonal matrix system. The population 
n is in a Boltzmann distribution at equilibrium, but not when out of equilibrium. Assuming the vibrational 

energy in the equilibrium form of a vibrational temperature, e v = e* v (T v ) , is inconsistent with the model! 

3) In the rate of change of the vibrational energy, a single time scale r is used in the model. This actually 
does not reflect the true collision processes. As we shall see later, the time scale for approaching 
equilibrium can be greatly different from the time scale leaving equilibrium. 

4) The model is invalid when dissociation occurs. 

Several models have been proposed to include the coupling of vibration and dissociation. The first coupled 
vibration-dissociation (CVD) model was proposed by Hammerling [6]. Heims [7] later derived more rigorous 
macroscopic equations using moments of the master equations. The CVD model was improved by Treanor and 
Marrone [8] and Marrone and Treanor [9] to allow the rate of vibrational excitation to be determined by the rate of 
dissociation. This coupled vibration-dissociation-vibration (CVDV) model was thus an attempt to account for a 
vibrational energy drain due to the dissociation process. Subsequently, Park developed an empirical two-temperature 
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model [10-14]. The model introduced an average temperature, which is a geometric average of the translational and 
vibrational temperatures, to dictate the dissociation rate and a correction factor to adjust the relaxation time at high 
temperatures for the Landau-Teller equation. Despite various attempts to improve the model’s physics, this type of 
vibrational modeling is only applicable for molecules, and does not apply to atoms or ions. It becomes essentially 
invalid at high temperature at which most molecules are dissociated into atoms. 

In looking at the aforementioned models, we see that all modeling had one thing in common - they were based 
on a single representation to cover the entire velocity or internal energy spaces. In contrast, modern numerical 
methods subdivide the physical space into elements or grids and employ piecewise representations directly in 
solving the governing equations, without further simplifications. In an effort to model various translational, thermal, 
and chemical non-equilibrium flow phenomena in a coupled and unified fashion, we are developing a new 
methodology [15] that blends some techniques used in the classical modeling [2-4, 16-17] and modem numerical 
methods [18,19]. 

Our modeling is based on the reconstmction of the state distribution function, which is expressed as a product 
of three space representation functions. Each representation is assumed to be a linear or nonlinear combination of 
some basis functions of its own space variables. The forms of the representations are based on maximum entropy, 
rational, and polynomial formulations, along with some other physical considerations. Each space is subdivided into 
multiple elements, groups, or bins, and piecewise representations are used in order to better describe the non- 
equilibrium gas mixture. The method of weighted residuals is applied to the microscopic governing equations in 
each space to form macroscopic moment equations. Unknowns are chosen to be either the moment integrals or their 
quadrature nodal values. The distribution function is reconstmcted in terms of the unknowns and then used to 
express the flux and collisional integral terms and the macroscopic rate coefficients in the moment equations. The 
modeling is completely physics-based, and its accuracy depends only on the assumed expression of the state 
distribution function and the number of subspaces used. In the simplest case, with all states being in one group, the 
current model recovers the classical Grad [16] and Levermore [17] models for translational non-equilibrium and 
produces a more general two-temperature model than the Landau-Teller [4] model for thermal and chemical non- 
equilibrium. In the other limiting case, in which each group contains only one state, the macroscopic moment 
equations reduce to the original microscopic governing equations. The current model also satisfies the translational- 
thermal-chemical equilibrium condition if there is detailed balancing and vice versa. 

In this paper, we present a study of our macroscopic modeling methodology applied to thermal and chemical 
non-equilibrium gases. Under these conditions, the microscopic governing equations are reduced to the master 
equations for this coupled system. This is one of few cases where we can practically solve the microscopic 
equations, even though this case is extremely stiff. Efficient numerical strategies for solving this system of nonlinear 
ordinary differential equations are discussed. We then derive the single-group and multi-group macroscopic 
equations and rate coefficients using our new methodology, with emphasis on the maximum entropy formulation. 
Finally, we show some numerical results and modeling validations. Macroscopic quantities such as density and 
internal energy from the solutions of the model are compared with the direct simulation of the microscopic 
equations. 


II. Microscopic Equations 


A. Master Equations 

In order to convey the basic idea of a new thermal and chemical non-equilibrium model, we consider a 
representative system involving nitrogen molecules undergoing internal excitation, de-excitation, dissociation, and 
recombination due to collisions with nitrogen atoms. The collisional processes can be expressed as 

Internal Excitation and De-excitation: N l 2 + N <=> N 2 + N (1) 

Dissociation and Recombination: N l 2 + N <=> N + N + N (2) 
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Here the superscript i or j is a generic index for an internal state of the molecules. If the energy of a molecular state i 
is above the dissociation energy, the molecule may undergo a non-collisional pre-dissociation process. In that case, 
iV asa collisional partner in Eq. (2) is not required. Note that the vibrational, rotational, and electronic motions of 
the nuclei and electrons of a realistic molecule are always coupled. Unlike classical modeling, we do not separate 
these modes in our modeling, even though each index i or j may include various quantum numbers associated with 
these motions. To simplify the notation, we list only a particular case here in which the internal state of atoms 
remains unchanged by the collision. Extensions to include all possible collisional and even radiative processes are 
straightforward, simply by summing all necessary terms over those processes. In this paper, the translational mode is 
assumed to be in equilibrium or to reach equilibrium instantaneously when the gas translational temperature is 
changed, satisfying the Maxwellian distribution in velocity. We further assume that the gas has no mean speed and 
that the state distribution function is not a function of physical position. The microscopic governing equations for 
the state population or number density are reduced to the master equations 

= ^(-k ,J n'n N +k J, n J n N )+(-k' d n'n N +k‘ r n N n N n N ) ■ co‘ , (3) 

= -2^V = - 2 ^(-K nin N + K n N n N n N ) ■ (4) 

dt i i 

Here n denotes the number density for the state i of N 2 and n N the number density for N . The symbols k ij and 
k jl represent the Maxwellian-distribution-based state-to-state excitation and de-excitation rate coefficients for the 
collisional process (1), and k l d and k l r the dissociation and recombination rate coefficients for the collisional process 
(2), respectively. These microscopic state-to-state rate coefficients satisfy the relations of detailed balancing. They 

are functions of the translational temperature T. We will also make use of the quantity !3 r = - — 

kT 


B. Direct Numerical Solver 


Equations (3) and (4) comprise a system of nonlinear ordinary differential equations representing the temporal 
evolution of the state populations of gas particles undergoing internal energy redistribution and chemical reaction 
processes to reach thermal and chemical equilibrium. The rate of change Jacobian can be obtained analytically 


A ik = ^- = -(W + £>A +k k ‘n N -2[V (-A;V +k i ‘n i )-k i d n i +3 k‘ r n N n N ] . 

r) n ^ ^ 


dri 


(5) 


In deriving Eq. (5), n N is not treated as an independent unknown. The algebraic version of Eq. (4), the conservation 
of element N, is used instead. The terms inside the brackets on the right hand side of Eq. (5) drop out if Eq. (4) is 
also solved. The Jacobian matrix is full and extremely stiff since it may consist of hundreds of thousands of internal 
states and the magnitude of n l may span over ten or even hundred orders of magnitude. The system therefore must 
be solved implicitly in time. The above equations can also be solve in a more efficient way by changing the 

unknowns to In , where g' denotes the degeneracy of the state i and n N . = V a* is the total molecular number 

n N 2 S i 

density. The governing equations for these new variables are then 


d n l 1 dn l _ 1 dn Ni 

dt n N g l n dt n N ^ dt 


co 


( 6 ) 


and the new Jacobian is 
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( 7 ) 


deb 1 


din 


n N,8 


n k n Ni da Y 
n N ^ - n k dn k 


The system has been studied previously [20, 21], using the ordinary differential equation (ODE) package 
LSODE [22]. In this paper we adopt a newer package SUNDIALS [23], in particular the ODE solver CVODE, for 
the time integration. The methods used in CVODE are variable-order, variable-step multistep methods. For stiff 
problems, CVODE employs the Backward Differentiation Formulas (BDFs), with order varying between 1 and 5. 
The nonlinear system is solved at each integration step using Newton iteration through either a direct (dense or 
banded) or a Krylov method linear solver. In the cases of the former, the Jacobian is fixed during the Newton 
iteration. When using the latter, the iteration is an Inexact Newton iteration, using the current Jacobian (through 
matrix-free products), in which the linear residual is nonzero but controlled. The SUNDIALS codes are written in C, 
but a set of FORTRAN interface routines are also supplied. The Jacobian matrix can be computed numerically 
within CVODE. However, using the analytical one can improve the CPU time by more than two times for this 
system. 


III. Macroscopic Modeling 


A. Local Representation and Reconstruction 

In a manner similar to many modern numerical methods [18,19] used in Computational Fluid Dynamics 
(CFD), our modeling is based on the reconstruction of the microscopic state number density n l , using some 

macroscopic quantities. Let e l denote the quantum energy level corresponding to the state i of N 2 . In our modeling, 
we assume n l can be expressed as a linear or nonlinear combination of some basis functions in terms of s l . As in 
CFD, the simplest choice is to use monomials or polynomials as the basis functions and express n l as a power 
series in e l 


n l = a + /3 s l + y (s l ) 2 + • • • . (8) 

The coefficients a , /3 , and y etc. are to be determined. This form of expression has a problem even at thermal 
equilibrium where n l follows the Boltzmann distribution. No matter how many terms we take in the power series, 
Eq. (8) simply cannot accurately represent an exponential function over the complete range of e l . One way to 
improve the accuracy is to use piecewise representations. We subdivide the internal energy space into multiple 
groups, similar to the grid cells or elements used in the physical space in CFD. However, the subdivision here is 
more flexible. It can be based on the values of e l , or vibrational and rotational quantum numbers. In each group g, 
we seek a representation 


n‘ =a g +p/ +y g (e‘f +■■■, iEg. (9) 

The accuracy and efficiency (smaller number of groups) can be further improved by choosing a better 
representation in each group. Based on the maximum entropy principle and the variational method [24], subject to 
the macroscopic constraints 

4 = S"' (£,)2 ’ (10) 

iEg iEg iEg 

we can derive the new representation 
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(11) 


In ?j- = a g +P g E‘ +y g (e‘) 2 +■■■, iEg. 

It is evident that Eq. (11) is a better representation than Eq. (9) at thermal equilibrium since the former can 
reproduce the Boltzmann distribution exactly using only one group with the constant and linear terms. The 
procedure to determine the coefficients, here in terms of the macroscopic constraints, is commonly referred as 
reconstruction in modern CFD methods. Substituting Eqs. (9) or (1 1) into (10), we can solve for the coefficients 

a g = a g( n g’ e g’fg’"'^ Pg = Pg( n g’ e g>J «’■■■)’ yg-yg( n g 9e g’fg’" m ')- (i2) 

Analytical solutions for Eq. (12) are not necessary. In fact, the form of the local representation and the macroscopic 
constraints can be arbitrary as long as the reconstruction is solvable. We would therefore prefer to choose ones that 
could make the solutions more accurate, and the evaluations of collisional terms in the macroscopic moment 
equations more efficient. This is an unexplored area that requires further investigation. 


B. Macroscopic Moment Equations and Rate Coefficients 

In order to reconstruct the microscopic state number density, we need the same number of macroscopic 
constraints as the number of terms used in the local representation. Using the method of weighted residuals, the 
equations governing the macroscopic constraints, or macroscopic moment equations, are obtained by multiplying the 
microscopic governing equations with appropriate weighting functions, substituting the reconstructed local 
representation for the microscopic state number densities, and then summing over all energy states in the group. 

We now give further details of the proposed model. We will employ the convention that the subscript g (or h ) 
is a dummy index of a group, while the superscript i (or j) is a dummy index of a state in the group g (or h). For 

simplicity we will omit the subscript g (or h) when referring to a particular state i (or j). Under this convention, n l 

and n J h are expressed as fl l and n J , k lJ gh and k£ as k lJ and k Jl . Symbols for all other variables follow the same 
convention. Equation (3) can now be expressed as 

< ^~ = y y (-£Vn v + k J, n J n N )+ (~k' d n'n N + k' r n N n N n N ). (13) 

^ h j&i 


Multiplying Eq. (13) by a weighting function ( e l ) m and summing over all microscopic states i in the group g , one 
obtains the m th moment equation of the master equation for the macroscopic constraint variables in Eq. (10) 


=J j '2'2{-k' J n'{e'T n N + k Jl n J {e'rn N 

i^gdt h i^gj^h 


) + 2 (-£>' (e‘ T n N + k' r (s' T n N n N n N ) . 

i^g 


(14) 


The right hand side of the above equation must be simplified to involve only macroscopic variables and rate 
coefficients. In this paper, we employ the maximum entropy formulation (1 1) to derive the generalized macroscopic 
rate equations. These equations can be expressed in a succinct manner by introducing moments of the partition 
function 


Q': = ^g'(eT 


+ Yg( £ ‘ ) 2+ -") 


i^g 


(15) 


The generalized macroscopic constraints or Eq. (10) can be replaced by the single equation 

«;=5>'Vr , 

i^g 


( 16 ) 
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where 


n„ = n„ 


<-f* 


(17) 


Note that the superscript m in the definitions of Q" 7 and n™ is not an exponent. Substituting Eq. (11) into (14), we 
obtain the relations 


and 


e ~ a s - Hl_ 

q: 


k gi 

< e! 


Thus, all macroscopic moment equations can be expressed by the single formulation 


dn 


u n — 

+ K n >N) + + K g n N n N n N ) . 


The generalized macroscopic rate coefficients are given by 


K >~, 




jE:h iE.g 
'dh i(Eg 

■ 

i^g 


-{P gE i +Yg {E i f + -) 


-(P h E j +Y h (E j f+-) 


- 0 6 g £ i +Y g {E i ) 2 + -) 


(18) 

(19) 

( 20 ) 

( 21 ) 

(22) 

(23) 

(24) 


In comparing Eqs. (20) with Eq. (3), we note that the macroscopic moment equations have the same form as the 
microscopic governing equations, retaining the original form of the collisional processes. They represent the time 
variation of the group population, internal energy, and other macroscopic quantities due to internal excitation and 
dissociation. The modeling involves only the interactions between the translational and internal modes. 

The modeling not only formulates the governing equations for macroscopic quantities, but also provides 
precise definitions of the macroscopic rate coefficients. Equations (21) - (24) indicate that these macroscopic 
coefficients are model-weighted sums of the corresponding microscopic rate coefficients. It is easy to show that 

( 25 ) 

but 


for Bil. 


(26) 


If the logarithm of the state distribution function is assumed to be linear in both the velocity and internal energy 
spaces, the macroscopic rate coefficients are then functions of two temperatures - the translational temperature and 


the internal temperature of the group. (The latter can be defined as — .) We note that the ratios of the macroscopic 

k/3g 


forward and backward rate coefficients are not equal to the (multi-temperature) equilibrium constant. This is 
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because multi-temperature “equilibrium” is not a true equilibrium state. This constant cannot be uniquely 
determined, but depends on which thermodynamic quantity is minimized or maximized in defining “equilibrium”. 


C. Special Cases for the Multi-Group Formulation 


The simplest way to divide internal energy states into groups is based on their energy levels. However, an 
optimum division into groups must consider all physical transition probabilities and may be guided by solving the 
complete master equations for a few cases. This will be carried out in a future study. Here we present a few special 
cases. 


a) One Energy State in a Group 

If any group contains only one microscopic energy state, then all higher order macroscopic moment equations 
for that group reduce to the group population equation (single energy level is a constant), and thus become 
redundant. If every group has only one energy state in its group, all macroscopic equations reduce to the original 
master equation; and all the macroscopic rate coefficients reduce to the microscopic state-to-state rate coefficients. 
This is limiting case provides a consistency check. Since the model must support this special case, it is not surprising 
that macroscopic group population and internal energy equations have the same form as the original master 
equations. 


b) Thermo-Chemical Equilibrium 

If the system is in thermo-chemical equilibrium, 

K = nf, p =($ „=& and r -y h - 0, 


(27) 


one can show that 


L^m m* _ j^m m* 
K gh n g ~ K hg H h 

js r m m* V m *' * 

K g d n g ~ K rg n N n N 


(28) 

(29) 


Here the superscript * denotes the translational-thermal-chemical equilibrium state. Equations (28)-(29) represent 
the detailed balancing for the groups. Conversely, if Eqs. (28)-(29) are satisfied, one can also show that the system is 
in thermo-chemical equilibrium. This is another limiting case that checks consistency. As we stated previously, one 
can show that, in general, 


and 


* K : g ^PrX 


ksp.ak 


*K™(l3 g ,l3 T )n: 


(30) 

(31) 


The equalities are also not true for the other group rate coefficients. One thus cannot relate the multi-group forward 
and backward rates to the equilibrium constant. 

c) All Energy States in One Group 

If we have only one group and all the energy states are in that group, we have /3 g = /3 h = /3 L * j3 r and n g = n h = n . 

Thus, we revert to a two-temperature model in terms of T and T , as we have defined 7] = — and T = — — . 

tp, kp r 

The first two terms on the right hand side of Eq. (20) for the group population cancel each other. This confirms the 
fact that the internal excitation and deexcitation do not change the total population. The zeroth moment of Eq. (20) 
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then reduces to the conventional chemical reaction equation involving dissociation and recombination. In the 
conventional use of the Two-Temperature model, the macroscopic rate constants are mostly based on the ground 
state rates, without considering excited state rates, and the temperature dependence is replaced with a geometric 
average of the translational and vibrational temperatures, T avg = T*T '~ S , where s is between 0.5 and 0.7 [3]. Here our 
model has a rigorous definition of the macroscopic rate coefficients based on all possible collisional processes. 

Unlike the total population equation, the first two terms on the right hand side of the first moment of Eq. (20) for 
the internal energy do not cancel each other since K l hg * K gh . This shows that the forward and backward energy 

transfer processes between the translational and internal modes have a different time scale. In contrast, the Landau- 
Teller equation is expressed in a relaxation form for which the two processes have the same time scale. Most models 
developed to improve the Landau-Teller equation were based on the assumption that the internal structure of 
molecules can be separated into macroscopic electronic, rotational, and vibrational modes. They thus include a very 
complicated modeling to deal with energy transfers among these modes and translational mode. Our model shows a 
much simpler form, involving only internal energy exchange with translational energy. More importantly, we make 
no assumptions at the microscopic level. All state-to-state collisional (and radiative) processes are allowed in the 
master equation. The model is applicable to both atoms and molecules and their ions. Since it is not justified to 
assign the vibrational temperature of a minor molecular species to the electronic temperature of a major atomic 
species, a simple extension of this case is to have separate internal temperatures for atomic and molecular species. 


d) Uniform Distribution Model 

If we set j3 = y g = • • • = 0 in each group, the current model reduces to the collisional-radiative (CR) model. 
Under this condition, the microscopic energy state populations in each group can be expressed as 

(32) 

i^g 



All macroscopic rate coefficients (21)-(24) have no dependence on ft . Since the group internal energy is not a 

constraint in this case, the formulation does not include any group internal energy equation. If one computes the 
group internal energy from other means, most likely this will be inconsistent with the second Eq. (10). In order for 
this special case to be a good approximation, one must have a large number of groups such that all energy levels in 
each group are very close together. Another drawback of this special case is that it does not reach thermal 
equilibrium. 


IV. Numerical Results and Modeling Validations 


Recently, scientists at NASA Ames Research Center have published a set of state-to-state collisional rate 
coefficients for nitrogen [25,26]. The dataset comprises 17 million collisional transitions between 9390 vibrational- 
rotational ( v,J) energy states of ground electronic state molecular nitrogen A 2 ('2p colliding with ground electronic 

state atomic nitrogen N( 4 S U ) . These include nonreactive and exchange-reactive vibrational and rotational 
excitations, dissociations, and noncollisional predissociation. The data set is used here to study the temporal 
evolution of the state populations of gas particles undergoing internal energy redistribution and chemical reaction 
processes to reach thermal and chemical equilibrium. Results of the study are then used to evaluate and validate our 
macroscopic models, in particular the prediction of macroscopic number density and internal energy. The state-to- 
state rate coefficients between nitrogen N 2 ( ] '2 + g ) molecules are currently being calculated. We will include 
molecule-molecule collisional processes in the future as those rates become available. 

Numerical calculations are carried out to study the non-equilibrium collisional process under isothermal 
conditions. Cold nitrogen molecules at room temperature 300 K, seeded with a small number of nitrogen atoms, are 
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suddenly heated, and the translational temperature of the gas is instantaneously brought up to a high temperature T . 
This drives the gas toward a very strong thermal and chemical non-equilibrium condition. The initial condition for 
pressure is set to 10,000 Pa, with a mixture of 95% N 2 and 5% N. Figure 1 shows the time evolution of the gas 
compositions (expressed as mole fractions) and internal energy of N 2 at T= 10,000 K. The figure indicates that 
dissociation takes place only after a period of incubation, approximately lqs. After the incubation period, it takes 
about 20 [is for the gas mixture to settle down and another 80 qs to reach chemical equilibrium. Figure 2 gives more 
details of the energy state population distributions at various time instants. The change of state populations with 
time shows how the translational energy is transferred to the internal energy through collisions. Initially, the state 
populations strongly depart from thermal equilibrium. The distribution of low-lying states is frozen at the initial 
temperature and the populations of high-lying states start to increase. At time t =10 ns, the high-lying state 
populations are nearly at equilibrium, although their distribution is still fairly broad. As relaxation proceeds, by the 
time t = 100 ns, the low-lying state populations are dissected into a few vibrational strands, while all other high- 
lying states have already reached equilibrium. This indicates that the energy states can be grouped for the purpose of 
macroscopic modeling. The entire thermal relaxation takes about ljis. The figure shows that the energy transfer 
process is almost complete before the onset of dissociation, indicating the two processes are independent at this 
temperature. This is made evident by the time evolution of the molecular internal energy shown in Figure 1. The plot 
reveals a plateau at time t = 1 |is for a few microseconds, indicating the gas is in thermal equilibrium. The internal 
energy then quickly decreases, showing the process of dissociation. At higher temperatures ( T > 10,000 - 
15,000K), dissociation occurs while the energy is being transferred. Figure 3 depicts the time evolution of the gas 
compositions (expressed as mole fractions) and internal energy of molecular nitrogen at T= 30,000 K. At this 
temperature, the gas is more energetic and therefore the dissociation occurs at an earlier time. The internal energy of 
N 2 never shows a plateau before it starts to decrease. This clearly indicates the two processes take place at the same 
time. Figure 4 shows the energy state population distributions at various time instants for T= 30,000 K. The 
variation of the distributions reaffirms that dissociation occurs while the energy is being transferred. Both Figs. (2) 
and (4) show the quasi-bound states (energy levels greater than the dissociation energy, 9.75eV for nitrogen) take a 
longer time than the bound states to reach thermal equilibrium. 

We next carry out the same numerical simulations using our macroscopic models. In this paper, we study only 
the simplest grouping strategy. The groups are defined by dividing the internal energy space uniformly. Figure 5 
depicts results for T = 10,000 K , using the uniform distribution model (/3 g = y g = ••• = 0 ). The figure compares the 

time evolution of the mole fraction of N atoms and the macroscopic internal energy of N 2 molecules for various 
groups. As shown in the figure, the model does not satisfy the thermal-chemical equilibrium condition, except for 
the limiting case of one state in each group. A converged solution is obtained if over 100 groups are used. However, 
these converged macroscopic solutions deviate slightly from the microscopic solution (labeled as Full Set). At the 
middle point, mole fraction equal to 0.5, the former is about 1.4 |is ahead of (faster) the latter. Using one group 
(each species is represented by one group), or even five or ten groups, the solution is far from the converged 
solution. We note that not all possible excitations and de-excitations were considered in the original data set [25]. In 
that data set, over 60% of the entries have no data. While the selection rule requires the jumps in rotational quantum 
number to be even for all inelastic collisions, exchange reactions do not have that selection rule, but have to be able 
to traverse the potential energy surface barrier. The averaging process in obtaining the macroscopic rate coefficients, 
Eqs. 21-24, simply reduces the number of zero entries. The percentage of zero entries in the macroscopic excitation 
and de-excitation rate coefficients is linearly proportional to the number of groups used, as shown in Figure 6. This 
changes the dynamics of the system, effectively making more groups participate in the collisional processes and 
therefore faster excitations and dissociations. Furthermore, the simple uniform internal energy grouping violates the 
selection rule. Many truly forbidden transitions thus become permissible. Figure 7 shows the same plots, but using a 
subset of data set involving only energy states whose vibrational quantum numbers v are less than or equal to 10. 
There are a total of 2845 internal energy states in this subset. The percentage of zero entries for the excitation rate 
coefficients is much smaller than the original one, and thus has less impact on the solution. From 100 groups and up, 
the macroscopic solutions match well to the microscopic solution. In order to further verify our argument regarding 
the influence of the zero entries to the converged solution, we filled the zero entries in the original microscopic rate 
coefficient dataset with the averaged values obtained from the 200-group macroscopic rate coefficients. The 
comparisons are shown in Figures 8a and 8b. While this is a purely mathematical exercise, the figures clearly show 
that the macroscopic solutions converge to the microscopic ones if over 100 groups are used for the uniform 
distribution model. 
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We have also included the results using the complete 9390 internal energy states and the original dataset for T 
= 30,000 K , as shown in Figure 9. The plots exhibit the same convergence pattern as the case for T = 10,000 K. 
However, the zero entries have a smaller impact in this case since the excitation and dissociation take place at the 
same time. We therefore observe a better agreement between the converged macroscopic solution and the 
microscopic solution. Due to the lack of a complete set of microscopic rate coefficients, we are unable to improve 
our model further using the uniform grouping. Nevertheless, the uniform distribution model shows that a converged 
solution can be obtained by solving approximately 100 macroscopic equations, a 99% reduction from 9390 
equations. In terms of CPU time, for the T= 10,000 K case using the analytical Jacobian, it takes about 18 hours to 
reach t = 10" 1 s on a MacBook Pro or a Linux workstation (single core) for the full set, but less than 1 second for the 
100 group case. 

We now discuss some results using the linear distribution model ( y = • • • = 0 ). Figures 10-11 depict the results 

for T = 10,000 K. The time evolution of the mole fraction for N atoms and the macroscopic internal energy of N 2 
molecules for various groups are given in Figure 10a and 10b, respectively. Comparing with Figure 5, the 
improvement of using the linear distribution model is evident. The plots show a much quicker convergence by 
increasing the number of groups. Even with one group or two groups, the macroscopic solutions are already very 
close to the exact microscopic ones. From three groups onward the macroscopic solutions are essentially the same. 
Again, the converged macroscopic solutions exhibit a slight deviation from the microscopic solution due to the lack 
of a complete set of microscopic excitation rate coefficients. Figure 11 shows the results for five groups, indicating 
the contribution of each group. Until about 100 ns, most of the contribution comes from group one, consisting of the 
low energy states. We have also carried out the same exercise using the filled microscopic rate coefficients for the 
linear distribution model. The results are given in Figure 12, showing that the macroscopic solutions converge to the 
microscopic solution even with only one group. This is because the internal excitation and dissociation are two 
independent processes for this temperature (T= 10,000 K). Figures 13 -14 show the same results as in Figures 10 - 
11 for a temperature of 30,000 K using the original dataset. Again, the plots exhibit the same convergence pattern as 
the case for T = 10,000 K. The macroscopic solutions also converge with three groups, but the converged solution 
has a better agreement with the microscopic solution since the incompleteness of the data set plays a smaller role in 
this case. Again, the CPU time for running the linear distribution model is less than 1 second. 


V. Concluding Remarks 


In this paper we presented a rigorous framework for modeling the thermal and chemical non-equilibrium gases. 
The model is based on the reconstruction of the state distribution function using the maximum entropy principle. 
The internal energy space is subdivided into multiple groups in order to better describe the non-equilibrium gases. 
The method of weighted residuals is applied to the microscopic equations to obtain macroscopic moment equations 
and rate coefficients. All the macroscopic moment equations and rate coefficients have the same forms as the 
microscopic counterparts. The modeling is completely physics-based, and its accuracy depends only on the assumed 
expression of the state distribution function and the number of groups used. The model makes no assumption at the 
microscopic level, and all possible collisional and radiative processes can be included. The model is applicable to 
both atoms and molecules and their ions. Several limiting cases are presented to show that the model recovers the 
classical two-temperature models if all states are in one group and the model reduces to the microscopic equations 
when there is only one state in each group. Numerical examples and model validations were carried out using the 
uniform and linear distribution models to study thermal and chemical non-equlibrium processes involving 
excitation, de-excitation, dissociation, and recombination in a nitrogen gas mixture. Results show that the original 
over nine thousand microscopic equations can be reduced to 2 macroscopic equations using 1 to 5 groups with 
excellent agreement. The computer time is decreased from 16 hours to almost none. In the future, we will investigate 
higher order distribution models. 

The modeling technique introduced here can be viewed as a new methodology applicable to a broad range of 
problems. We envision this can be applied to translational non-equilibrium flows as well. We are carrying out a 
study to model various translational, thermal, and chemical (and possibly radiative) non-equilibrium flow 
phenomena in a coupled and unified fashion. As a final remark, we would like to re-emphasize that the key quantity 
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for much of high-temperature gas physics and radiation modeling is the state-to-state transition probability or cross 
section. It simply makes no sense to carry out any numerical simulations if one could not obtain the gas properties 
and rate coefficients correctly. Currently, many state-of-the-art CFD design tools for space vehicles still depend on 
tuning parameters or ad hoc fixes in some areas. We do have the capability to compute those cross sections in 
computational chemistry. As we are entering a new era of space flight that the speed is beyond lunar return, now is 
the time to start computing those cross sections and carrying out modeling in a more physics-based fashion! 
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Figure 1. Mole fractions of N and N 2 and internal energy of N 2 at T = 10,000K 
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Figure 2. Energy state population distributions at various time instants, T= 10,000 K. 
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Figure 3. Mole fractions of N and N 2 and internal energy of N 2 at T = 30,000K 
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Figure 4. Energy state population distributions at various time instants, T=3 0,000 K. 
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Figure 5. Comparison of time evolution for T= 10,000 K using uniform distribution model 
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Figure 6. Percentage of zero entries in excitation rate coefficients. 
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Figure 7. Comparison of time evolution for v < 10 using uniform distribution model 
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Figure 8. Comparison of time evolution using uniform distribution model with filled data set. 
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Figure 9. Comparison of time evolution for T= 30,000 K using uniform distribution model 
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Figure 10. Comparison of time evolution for T= 10,000 K using linear distribution model 
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Figure 11. Individual group quantity T= 10,000 K using linear distribution model. 
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Figure 12. Comparison of time evolution using uniform distribution model with filled data set. 
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Figure 10. Comparison of time evolution for T= 30,000 K using linear distribution model 
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Figure 11. Individual group quantity T= 10,000 K using linear distribution model. 
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